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Abstract — Existing approaches to image reconstruction in 
pliotoacoustic computed tomography (PACT) with acoustically 
heterogeneous media are limited to weakly varying media, are 
computationally burdensome, and/or cannot effectively mitigate 
the effects of measurement data incompleteness and noise. In 
this work, we develop and investigate a discrete imaging model 
for PACT that is based on the exact photoacoustic (PA) wave 
equation and facilitates the circumvention of these limitations. A 
key contribution of the work is the establishment of a procedure 
to implement a matched forward and backprojection operator 
pair associated with the discrete imaging model, which permits 
application of a wide-range of modern image reconstruction 
algorithms that can mitigate the effects of data incompleteness 
and noise. The forward and backprojection operators are based 
on the k-space pseudospectral method for computing numerical 
solutions to the PA wave equation in the time domain. The 
developed reconstruction methodology is investigated by use of 
both computer-simulated and experimental PACT measurement 
data. 

Index Terms — Photoacoustic tomography, optoacoustic tomog- 
raphy, thermoacoustic tomography, 
iterative image reconstruction, acoustic heterogeneity 



I. Introduction 

Photoacoustic computed tomography (PACT), also known 
as optoacoustic or thermoacoustic tomography, is a rapidly 
emerging hybrid imaging modality that combines optical 
image contrast with ultrasound detection. [l]-[4] In PACT, 
the to-be-imaged object is illuminated with a pulsed optical 
wavefield. Under conditions of thermal confinement [2], [5], 
the absorption of the optical energy results in the generation of 
acoustic wavefields via the thermoacoustic effect. These wave- 
fields propagate out of the object and are measured by use of 
wide-band ultrasonic transducers. From these measurements, a 
tomographic reconstruction algorithm is employed to obtain an 
image that depicts the spatially variant absorbed optical energy 
density distribution within the object, which will be denoted by 
the function A{r). Because the optical absorption properties 
of tissue are highly related to its hemoglobin concentration 
and molecular constitution, PACT holds great potential for a 
wide-range of anatomical, functional, and molecular imaging 
tasks in preclinical and clinical medicine [3], [6]-[9]. 
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The majority of currently available PACT reconstruction 
algorithms are based on idealized imaging models that assume 
a lossless and acoustically homogeneous medium. However, 
in many applications of PACT these assumptions are violated 
and the induced photoacoustic (PA) wavefields are scattered 
and absorbed as they propagate to the receiving transducers. 
In small animal imaging applications of PACT, for example, 
the presence of bone and/or gas pockets can strongly perturb 
the photoacoustic wavefield. Another example is transcranial 
PACT brain imaging of primates [10], in which the PA wave- 
fields can be strongly aberrated and attenuated [11]-[13] by the 
skull. In these and other biomedical applications of PACT, the 
reconstructed images can contain significant distortions and 
artifacts if the inhomogeneous acoustic properties of the object 
are not accounted for in the reconstruction algorithm. 

Several image reconstruction methods have been proposed 
to compensate for weak variations in a medium's speed-of- 
sound (SOS) distribution [14]-[16]. These methods are based 
on geometrical acoustic approximations to the PA wave equa- 
tion, which stipulate that the PA wavefields propagate along 
well-defined rays. For these ray-based propagation models to 
be valid, variations in the SOS distribution must occur on 
length scales that are large compared to the effective acoustic 
wavelength. These assumptions can be violated in preclinical 
and clinical applications of PACT. To compensate for strong 
SOS variations, a statistical approach has been proposed [17] 
to mitigate the artifacts in the reconstructed images caused 
by the wavefront distortions by use of a priori information 
regarding the acoustic heterogeneities. However, this method 
neglected variations in the medium's mass density and the 
effects of acoustic attenuation. 

A few works have reported the development of full-wave 
PACT reconstruction algorithms that are based on solutions to 
the exact PA wave equation [18]-[23]. While these methods 
are grounded in accurate models of the imaging physics 
and therefore have a broader domain of applicability than 
ray-based methods, they also possess certain practical lim- 
itations. Finite element methods (FEMs) have been applied 
for inverting the PA wave equation in both the time and 
temporal frequency domains [18], [19]. However, a very large 
computational burden accompanies these methods, which is 
especially problematic for three-dimensional (3D) applications 
of PACT. Image reconstruction methods based on time-reversal 
(TR) are mathematically exact in their continuous forms 
in homogeneous media for the 3D case [20]. While these 
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methods possess significantly lower computational burdens 
then FEM-based approaches, they possess other limitations 
for use with practical PACT applications. For example, TR 
methods are predicated upon the assumption that the measured 
PA signals are densely sampled on a measurement surface that 
encloses the object, which is seldom achievable in biomedical 
applications of PACT. More recently, a Neumann series-based 
reconstruction method has been reported [22], [23] for media 
containing SOS variations that is based on a discretization of a 
mathematically exact inversion formula. The robustness of the 
method to practical sparse sampling of PA signals, however, 
has not been established. 

In this work, we develop and investigate a full-wave ap- 
proach to iterative image reconstruction in PACT with media 
possessing inhomogeneous SOS and mass density distributions 
as well as acoustic attenuation described by a frequency 
power law. The primary contributions of the work are the 
establishment of a discrete imaging model that is based on 
the exact PA wave equation and a procedure to implement 
an associated matched discrete forward and backprojection 
operator pair The availability of efficient numerical procedures 
to implement these operators permits a variety of modern itera- 
tive reconstruction methods to be employed that can effectively 
mitigate image artifacts due to data incompleteness, noise, 
finite sampling , and modeling errors. Specifically, the k-space 
pseudospectral method is adopted [21] for implementing the 
forward operator and a numerical procedure for implementing 
the exact adjoint of this operator is provided. The k-space 
pseudospectral method possesses significant computational 
advantages over real space finite-difference and finite-element 
methods, as it allows fewer mesh points per wavelength and 
allows larger time steps without reducing accuracy or intro- 
ducing instability [24]. An iterative image reconstruction algo- 
rithm that seeks to minimize a total variation (TV)-regularized 
penalized least squares (PLS) cost function is implemented 
by use of the developed projection operators and investigated 
in computer-simulation and experimental studies of PACT in 
inhomogeneous acoustic media. Also, the performance of this 
algorithm is compared to that of an existing TR method. 

The paper is organized as follows. In Section II, the 
salient imaging physics and image reconstruction principles 
are briefly reviewed. The explicit formulation of the discrete 
imaging model is described in Section III. Section IV gives a 
description of the numerical and experimental studies, which 
includes the implementation of the forward and backprojection 
operators, and the iterative reconstruction algorithm. The nu- 
merical and experimental results are given in Section V. The 
paper concludes with a summary and discussion in Section 
VI. 

II. Background 

Below we review descriptions of photoacoustic wavefield 
generation and propagation in their continuous and discrete 
forms. The discrete description is based on the k-space pseu- 
dospectral method [21], [24], [25]. We present the pseu- 
dospectral k-space method by use of matrix notation, which 
facilitates the establishment of a discrete PACT imaging model 



in Section III. We also summarize a discrete formulation of 
the image reconstruction problem for PACT in acoustically 
inhomogeneous media. Unless otherwise indicated, lowercase 
and uppercase symbols in bold font will denote vectors and 
matrices, respectively. 

A. Photoacoustic wavefield propagation: Continuous formu- 
lation 

Let p(r, i) denote the thermoacoustically-induced pressure 
wavefield at location r e M'^ and time t > {). Additionally, 
let A{y) denote the absorbed optical energy density within the 
object, r(r) denote the dimensionless Grueneisen parameter, 
u(r,i) = (w^(r, i), w^(r, t), u^(r, t)) denote the vector- valued 
acoustic particle velocity, co(r) denote the medium's SOS 
distribution, and p{r,t) and po{r) denote the distributions of 
the medium's acoustic and ambient densities, respectively. The 
object function A{r) and all quantities that describe properties 
of the medium are assumed to be represented by bounded 
functions possessing compact supports. 

In many applications, acoustic absorption is not negligible 
[21], [26]-[29]. For a wide variety of lossy materials, including 
biological tissues, the acoustic attenuation coefficient a can be 
described by a frequency power law of the form [30] 

a{vJ) = ao{r)fy, (1) 

where / is the temporal frequency in MHz, ao is the 
frequency-independent attenuation coefficient in dB MHz^^ 
cm~^, and y is the power law exponent which is typically in 
the range of 0.9-2.0 in tissues [31]. 

In a heterogeneous lossy fluid medium in which the acoustic 
absorption is described by the frequency power law, the 
propagation of p{r, t) can be modeled by the following three 
coupled equations [21], [32] 

d 1 

-u(r,t)- r^Vp(r,t), (2) 

ot pq[y) 

^p(r,i) = -po(r)V-u(r,<), (3) 

p(r,t)=co(rf{l-Mr)|(-V2)''/'-^ 

-,7(r)(-V2)(^'-i)/2}p(r,t), 

subject to the initial conditions: 

Po(r)=p(r,i)|t=o-r(r)A(r), u(r, t)|t=o = 0. (5) 

where the quantities /i(r) and 7](r) describe the acoustic 
absorption and dispersion proportionality coefficients that are 
defined as 

/x(r) = -2aoCo(r)^-\ 7;(r) = 2aoCo(r)^tan(7ry/2). 

(6) 

Note that acoustic absorption and dispersion are modeled by 
the second and third terms in the bracket, which employ two 
lossy derivative operators based on the fractional Laplacian to 
separately account for the acoustic absorption and dispersion 
in a way that is consistent with Eqn. (1). When acoustic 
attenuation can be neglected, /i(r) =0 and 7;(r) = 0, and 
Eqn. (4) reduces to 

p(r,t) =co(r)V(r,t). (7) 
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B. Photoacoustic wavefield propagation: Discrete formulation 

The k-space pseudospectral method can be employed to 
propagate a photoacoustic wavefield forward in space and time 
by computing numerical solutions to the coupled equations 
described by Eqn. (2), (3), (4), and (5). This method can be sig- 
nificantly more computationally efficient than real space finite- 
element and finite-difference methods because it employs the 
fast Fourier transform (FFT) algorithm to compute the spatial 
partial derivatives and possesses less restrictive spatial and 
temporal sampling requirements. Applications of the k-space 
pseudospectral method in studies of PACT can be found in 
references [10], [13], [21], [24]. 

The salient features of the k-space pseudospectral method 
that will underlie the discrete PACT imaging model are 
described below. Additional details regarding the application 
of this method to PACT have been published by Treeby and 
Cox in references [21], [24]. Let ri,--- ,rAr e R"^ specify 
the locations of the N = N1N2N3 vertices of a 3D Cartesian 
grid, where Ni denotes the number of vertices along the i- 
th dimension. Additionally, let niAt, m e 7L* , At e M+, 
denote discretized values of the temporal coordinate t, where 
Z* and M+ denote the sets of non-negative integers and 
positive real numbers. The sampled values of p{r, t = mAt) 
and ■u*(r, t = mAt), i = 1,2 or 3, corresponding to spatial 
locations on the 3D Cartesian grid will be described by the 
3D matrices P™ and U^, respectively, where the subscript m 
indicates that these quantities depend on the temporal sample 
index. Unless otherwise indicated, the dimensions of all 3D 
matrices will be Ni x N2 x N3. Lexicographically ordered 
vector representations of these matrices will be denoted as 



and 



= (u'(ri,77iAt),--- ,u'{rN,mAt)f, (8) 



ip{ri,mAt),--- ,p{rN,mAt)f. (9) 



The sampled values of the ambient density po{r) and squared 
SOS distribution CQ(r) will be represented as 



Q = diag(po(ri 



,Po(rjv)), 



and 



diag(cg(ri),--- ,cl{rN)), 



(10) 



(11) 



where diag(ai, a^?) defines a diagonal 2D matrix whose 
diagonal entries starting in the upper left corner are ai, ajy. 

In the k-space pseudospectral method, the ID discrete 
spatial derivatives of the sampled fields with respect to the 
i-th dimension (i = 1,2, or 3) are computed in the Fourier 
domain as 



VMatp 



F-i{iK'oKoF{P„}}, 



and 



vru:„^F-io-K'oKoF{u:„}}, 



(12) 



(13) 



where j = ^/— T, the superscript 'Mat' indicates that the 
ID discrete derivative operator V^^' acts on a 3D matrix, F 
and F^^ denote the 3D forward and inverse discrete Fourier 



transforms (DFTs), and o denotes Hadamard product. The 
elements of the 3D matrix K* (i = 1, 2, 3) are given by 

ni — 1 



(14) 



where = 1, • • • ^Ni (i = 1, 2, 3), and Li denotes the length 
of the spatial grid in the i-th dimension. 

The 3D matrix k = sinc(i AtCminK) is the k-space operator. 



where sinc(x) = £i!l(£i 
3D matrix defined as 



, Cmin is the minimum of co(r), K is a 



K 



\ 1=1 



(15) 



and the sine function and square root function are both 
element- wise operations. 

Consider the operators €>f ^' and 'i'^^' that are defined as 



*Matp 

^ A ^ r. 



-AfQ-iVf^'P„ 



and 



(16) 



(17) 



It will prove convenient to introduce the NxN matrices $i 
and that act on the vector representations of the matrices 
Pm and U^j, respectively. Specifically, $i and are defined 
such that ^iPm and ^iUj^ are lexicographically ordered 
vector representations of the matrices $f ^'P^ and \I'^^'UJ„, 
respectively. In terms of these quantities, the discretized forms 
of Eqn. (2), (3), and (4) can be expressed as 



Pm+1 = Pm + 



(18) 
(19) 



where pl^ is an x 1 vector whose elements are defined to 
be zero for m = 0, and 

3 

p,n+l = C Y,{p'm+1 + + ^P'm+l}- (20) 

1=1 

The quantities Au^^^ ™d ^Pm+i ^1"- (20) represent the 
absorption and dispersion terms in the equation of state. They 
are defined as lexicographically ordered vector representations 
of A^'"U^j_|_i and B^'"N^_|_i, which are defined in analogy 
to Eqn. (4) as 



1=1 



(21) 



B^'-'Nj^+i ^ TjF-i |k«-1f{ ^N:„+i}| , (22) 

where N^^j^ is the 3D matrix form of pl^, and p, and r] are 
defined as 



p = diag(/io(ri), • • • ,^o(rjv)), 
77 = diag(?/o(ri), ■ ■ ■ ,?7o(rjv)), 



(23) 
(24) 
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and ^ and ^ are powers of K that are computed on 
an element-wise basis. 

C. The image reconstruction problem 

Here, for simplicity, we neglect the acousto-electrical im- 
pulse response (EIR) of the ultrasonic transducers and assume 
each transducer is point-like. However, a description of how to 
incorporate the transducer responses in the developed imaging 
model is provided in Appendix-A. With these assumptions, 
we can define p„j = (p(rf , mAi), • • • ,p(r^, mAt))"'" as 
the measured pressure wavefield data at time t ~ mAt 
(to = 0, • • • , M — 1), where M is the total number of time 
steps and e (l = 1 , • • • , L) denotes the positions of 
the L ultrasonic transducers that reside outside the support 
of the object. The PACT image reconstruction problem we 
address is to obtain an estimate of po{r) or, equivalently, 
A{r), from knowledge of p™, to = O,--- ,A/ — 1, co(r), 
Pa{r), cto{r), and y. The development of image reconstruction 
methods for addressing this problem is an active area of 
research [10], [20], [22], [24], [33]. Note that the acoustic 
parameters of the medium can be estimated by use of adjunct 
ultrasound tomography image data [34]-[36] and are assumed 
to be known. The effects of errors in these quantities on the 
accuracy of the reconstructed PACT image will be investigated 
in Section IV. 

The discrete form of the imaging model for PACT can be 
expressed generally as 



P = Hpo, 



where the LM x 1 vector 



Po 
Pi 

Pm-1 



(25) 



(26) 



represents the measured pressure data corresponding to all 
transducer locations and temporal samples, and the x 1 
vector Po is the discrete representation of the sought after 
initial pressure distribution within the object (i.e., Eqn. (9) 
with m = 0). The LAI x N matrix H represents the discrete 
imaging operator, also referred to as the system matrix. 

The image reconstruction task is to determine an estimate 
of Po from knowledge of the measured data p. This can 
be accomplished by computing an appropriately regularized 
inversion of Eqn. (25). When iterative methods are employed 
to achieve this by minimizing a penalized least squares cost 
function [37], the action of the operators H and its adjoint 
must be computed. Methods for implementing these operators 
are described below. 

III. Explicit formulation of discrete imaging 

MODEL 

The k-space pseudospectral method for numerically solving 
the photoacoustic wave equation described in Section II-B will 
be employed to implement the action of the system matrix H. 
In this section, we provide an explicit matrix representation of 
H that will subsequently be employed to determine H^. 



Equations (18) - (20) can be described by a single matrix 
equation to determine the updated wavefield variables after a 
time step At as 

v„i+i = Wv,„, (27) 



where v„ 



1 Ll™, 



^m' Pm 1 Pm 



Pm^PniV is a 7A^ X 1 
vector containing all the wavefield variables at the time step 
771 At. The 7N x IN propagator matrix W is defined as 



W = 

OnxN 
OnxN 

*1 
OnxN 
OnxN 

Di 



OnxN 

OnxN 

OnxN 
T>2 



OnxN 
OnxN 
^NxN 
OnxN 
OnxN 
*3 



OnxN 
OnxN 
OnxN 

OnxN 
OnxN 

E 



OnxN 
OnxN 
OnxN 
OnxN 
'i-NxN 
OnxN 

E 



OnxN 
OnxN 
OnxN 
OnxN 
OnxN 
^NxN 

E 



*1 

G 

(28)' 



where ~ C(A + ^^ + B*,) (i = 1, 2, 3), E = C + CB, 

3 

G = CY, A#, + (I + B)*i*i, Inxn is the iVx identity 

1=1 

matrix, and O^vxAf is the N x N zero matrix. Recall that S&i 
was defined below Eqn. (17). 

The wavefield quantities can be propagated forward in time 
from t ^ to t = (M ~ l)At as 



vo 

Vl 

VM-1 



■ ■ 



•Ti 



Vo 
OtNxI 







7JVxl 



where the 7NM x 7A^M matrices {m = l,- 
are defined in terms of W as 



(29) 



, M - 1) 



f-lNxTN 



'TNxTN 







07Arx7Af ••• ItJVxTJV 
QlNxlN ■■■ W 
0(M-m-l)-7Afxm-77V 



(m+l)-7Nx(M~m)-7N 



>{M-m-l)-7Nx(M-m)-7N 

(30) 



with W residing between the {7N{m — 1) + l)-th to 7iVTO-th 
rows and the (7A^to + l)-th to 7A^(to + l)-th columns of T^. 

From the equation of state in Eqn. (7) and initial conditions 
Eqn. (5), the vector (vo,0, • • • ,0)^ can be computed from 
the initial pressure distribution po as 



Vo 
O7NXI 



J7Nxl 



ToPo, 



where 



To= (x,07 



TNxN, ■ ■ ■ , OjNxn) , 



(31) 



(32) 



T = (OatxAT, O^fxTV, OtvxAt, o*-" ^'o^ ^'Q^ ^i'^Nxn)'^, 

^ ^ ^ (33) 
and Po is the initial pressure distribution as defined by Eqn. 
(9) with TO = 0. 
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In general, the transducer locations at which the PA 
data p are recorded will not coincide with the vertices of 
the Cartesian grid at which the values of the propagated field 
quantities are computed. The measured PA data p can be 
related to the computed field quantities via an interpolation 
operation as 

vo 



p = S 



Vl 



VA/-1 



(34) 



where 



S = 







lLx7N 



>Lx7N 



Lx7N ■ ■ 




>Lx7N 



Olx7N 
0lx7N 





(35) 



Here, = [si,--- ,8^]"^, where s; (/ — I,-- - ,i) is a 
1 X 7N row vector in which all elements are zeros except 
the 4 corresponding to acoustic pressure at 4 grid nodes 
r;,i, i"i,2, ri,3, r;,4 that are nearest to the transducer location rf. 
In other words, these 4 entries are interpolation coefficients 
to compute the acoustic pressure at the Z-th transducer, and 
their values are given by the barycentric coordinates of rf 
with respect to rj.i, r; 2, ri,3, r;.4, which are determined by 
Delaunay triangulation [38]. 

By use of Eqns. (29), (31), and (34), one obtains 



p = ST 



■TiToPo. 



(36) 



Finally, upon comparison of this result to Eqn. (25), the 

sought-after explicit form of the system matrix is identified 
as 

H = STm-i---TiTo. (37) 

Commonly employed iterative image reconstruction meth- 
ods involve use of a backprojection matrix that corresponds 
to the adjoint of the system matrix. Since H contains real- 
valued elements in our case, is equivalent to the transpose 
H^. According to Eqn. (37), the explicit form of is given 
by 



tJtt...tt_iST. 



(38) 



The implementations of H and are described in Section 
IV- A. Note that, although the descriptions of H and above 
are based on the 3D PA wave equation, the two-dimensional 
formulation is contained as a special case. 



A. Implementation of the forward and backprojection opera- 
tors 

The k-space pseudospectral method for numerically solving 
the photoacoustic wave equation has been implemented in the 
MATLAB k-Wave toolbox [39]. This toolbox was employed 
to compute the action of H. To prevent acoustic waves from 
leaving one side of the grid and re-entering on the opposite 
side, an anisotropic absorbing boundary condition called a 
perfectly matched layer (PML) was employed to enclose 
the computational grids. The performance of the PML was 
dependent on both the size and attenuation of the layer A PML 
thickness of 10 grid points, together with a PML absorption 
coefficient of 2 nepers per meter, were found to be sufficient 
to reduce boundary reflection and transmission for normally 
incident waves [40], [41] and were employed in this study. To 
accurately and stably model wave propagation, the temporal 
and spatial steps were related by the Courant-Friedrichs-Lewy 
(CFL) number as [25], [39] 



At < 



CFLAr„ 



(39) 



where the Ar^in is the minimum grid spacing, and a CFL 
number of 0.3 typically provides a good compromise between 
computation accuracy and speed [39], [40]. A more detailed 
description of the implementation of the k-space pseudospec- 
tral method can be found in Refs. [39], [40]. 

The action of the backprojection matrix on the measured 
pressure data p was implemented according to Eqn. (38). It 
can be verified that p'^P = H^p can be computed as 



0'^Pa/-i, 



Af - I,-- 



.bp 



(40) 
,1 (41) 
(42) 



Since and r are both sparse matrices that can be stored 
and transposed, 0^p„i and r^v^ can be readily computed. 
Most of block matrices in the propagator matrix W are zero 
or identity matrices. Therefore, to compute W^v™, we only 
need to compute the actions of transposed non-trivial block 
matrices in W. To incorporate the PML boundary condition, 
both W and should be modified as described in Ref. 
[40]. 



IV. Descriptions of numerical and experimental 

STUDIES 

Numerical studies were conducted to demonstrate the ef- 
fectiveness and robustness of the proposed discrete imaging 
model in studies of iterative image reconstruction from incom- 
plete data sets in 2D and 3D PACT. Specifically, the system 
matrix and its adjoint, as formulated in Section III, were 
employed with an iterative image reconstruction algorithm that 
was designed to minimize a PLS cost function that contained 
a total variation (TV) penalty term. The performance of the 
reconstruction algorithm was compared to an existing TR- 
based reconstruction algorithm. 



B. Reconstruction algorithms 

By use of the proposed discrete imaging model and methods 
for implementing H and H^, a wide variety of iterative image 
reconstruction algorithms can be employed for determining 
estimates of po. In this work, we utilized an algorithm that 
sought solutions of the optimization problem 



Po 



arg mm 

po>0 



p-Hpoll' +A|po|tv, 



(43) 



where A is the regularization parameter, and a non-negativity 
constraint was employed. For the 3D case, the TV-norm is 
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defined as 

N 

|Po|TV = ^{([Po]n-[Po]„-)' + 

n=l 

([P0]n-[P0]„-)' + ([P0]„-[P0]„-f}^ 

(44) 

where [po]n denotes the n-th grid node, and [po]„-, 
[Po]„- J [Po]„- ^6 neighboring nodes before the n-th node 
along the first, second and third dimension, respectively. The 
fast iterative shrinkage/thresholding algorithm (FISTA) [42], 
[43] was employed to solve Eqn. (43), and its implementation 
is given in Appendix-B. The regularization parameter A was 
empirically selected to have a value of 0.001 and was fixed 
for all studies. 

A TR image reconstruction algorithm based on the k-space 
pseudospectral [21] method was also utilized in the studies 
described below. The TR reconstruction algorithm solves the 
discretized acoustic Eqns. (18) - (20) backward in time subject 
to initial and boundary conditions as described in reference 
[21]. The parameters of the PML boundary condition were 
the same with the ones employed in our system matrix 
construction. 

For both algorithms, images were reconstructed on a uni- 
form grid of 512 x 512 pixels with a pitch of 0.2 mm for the 
2D simulation studies and on a 256 x 256 x 128 grid with 
a pitch of 0.4 mm for the 3D studies. All simulations were 
computed in the MATLAB environment on a workstation that 
contained dual hexa-core Intel(R) Xeon(R) E5645 CPUs and 
a NVIDIA Tesla C2075 GPU. The GPU was equiped with 
448 1.15 GHz CUDA Cores and 5 GB global memory. The 
Jacket toolbox [44] was employed to perform the computation 
of Eqns. (18) - (20) and (40) - (42) on the GPU. 

C. Computer-simulation studies of 2D PACT 

Scanning geometries: Three different 2D scanning geome- 
tries were considered to investigate the robustness of the 
reconstruction methods to different types and degrees of data 
incompleteness. A 'full-view' scanning geometry utilized 180 
transducers that were evenly distributed on a circle of radius 40 
mm. A 'few-view' scanning geometry utiUzed 60 transducers 



that were equally distributed on the circle. Finally, a 'limited- 
view' scanning geometry utilized 90 transducers that were 
evenly located on a semi-circle of radius 40 mm. 

Numerical phantoms: The two numerical phantoms shown 
in Fig. l-(a) and (b) were chosen to represent the initial 
pressure distributions po in the 2D computer-simulation stud- 
ies. The blood vessel phantom shown in Fig. l-(a) was 
employed to investigate the robustness of the reconstruction 
methods with respect to different types and degrees of data 
incompleteness mentioned above. The low contrast disc phan- 
tom displayed in Fig. l-(b) was employed to investigate the 
robustness of the reconstruction methods with respect to errors 
in the SOS and density maps introduced below. 

Measurement data: Assuming ideal point-like transducer 
and neglecting the transducer EIR and acoustic attenuation, 
simulated pressure data corresponding to the numerical phan- 
toms were computed at the transducer locations by use of 
the k-space pseudospectral method for the 3 measurement 
geometries. To avoid committing an 'inverse crime' [45], a 
1024 X 1024 grid with a pitch of 0.1 mm was employed 
in this computation. A total of 20,000 temporal samples 
were computed at each transducer location with time step 
At = 30 ns, all of which were employed by the TR image 
reconstruction method. However, only the first 1,500 temporal 
samples were employed by the iterative reconstruction method. 
The same procedure was repeated for noisy pressure data, 
where 3% (with respect to maximum value of noiseless data) 
additive white Gaussian noise (AWGN) was added to the 
simulated pressure data. 

Investigation of systematic errors: The SOS and density 
maps employed in the simulation studies were representative 
of a monkey skull [10]. The dimensions of the skull were 
approximately 7 cm x 6 cm, and its thickness ranges from 2 
to 4 mm. Figure 2-(a) and (b) show a transverse slice of the 
SOS and density maps, which were used in the 2D simulations. 

Since errors in the estimated SOS and density maps are 
inevitable regardless in how they are determined, we inves- 
tigated the robustness of the reconstruction methods with 
respect to the SOS and density map errors, which were 
generated in two steps. First, 1.3% (with respect to maximum 
value) uncorrected Gaussian noise with mean value of 1.7% 
of the maximum value was added to the SOS and density 
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Fig. 2. A slice of tlie SOS (a) and density (b) map deduced from the X-ray 
CT data of a monkey skull. Panel (c) and (d) display profiles of the SOS and 
density maps along the 'X'-axis indicated in Fig. 2, respectively. Red dashed 
lines ai'e the profiles of the assumed maps, whereas the blue solid lines are 
the profiles of maps with errors. 

maps to simulate inaccuracy of the SOS and density values. 
Subsequently, the maps were shifted by 7 pixels (1.4 mm) to 
simulate a registration error Figure 2-(c) and (d) show profiles 
of the SOS and density maps with those errors along the 'X'- 
axis indicated by the arrows in Fig. 2-(a) and (b), respectively. 

D. Computer-simulation studies of 3D PACT 

Because PACT is inherently a 3D method, we also con- 
ducted 3D simulation studies to evaluate and compare the 
iterative reconstruction method and the TR method. As in 
the 2D studies described above, the 3D SOS and density 
maps were representative of a monkey skull. A 3D blood 
vessel phantom was positioned underneath the skull to mimic 
the blood vessels on the cortex surface. To demonstrate this 
configuration. Figure l-(c) shows the overlapped images of the 
3D phantom and the skull. The assumed scanning geometry 
was a hemispherical cap with radius of 46 mm, and 484 
transducers were evenly distributed on the hemispherical cap 
by use of the golden section spiral method [46]. The pressure 
data were computed on a 512 x 512 x 256 grid with a pitch 
of 0.2 mm and a time step At — 30 ns. The simulated 
pressure data were then contaminated with 3% AWGN. The 
TR reconstruction method employed 2,000 temporal samples 
at each transducer location, whereas the iterative method 
employed 1,000 samples. 

E. Studies utilizing experimental data 

Since the acoustic absorption and dispersion were modeled 
by the system matrix, the iterative method can naturally 
compensate for absorption and dispersion effects during recon- 
struction. To demonstrate the compensation for those effects, 
images were reconstructed by use of the iterative method with 




Fig. 3. A photograph of the pencil leads held in agar and sun'ounded by an 
acrylic cylindrical shell. 



experimental data obtained from a well-characterized phantom 
object that is displayed in Fig. 3. The phantom contained 6 op- 
tically absorbing structures (pencil leads with diameter 1 mm) 
embedded in agar. These structures were surrounded by an 
acrylic cylinder, which represents the acoustic heterogeneities 
and absorption in the experiments. The cylinder had inner and 
outer radii of 7.1 and 7.6 cm, respectively, and a height of 3 
cm. The density and SOS of the acrylic were measured and 
found to be 1200 kg m"'^ and 3100 m s~^, and the estimated 
acoustic absorption parameters were found to be ao = 1-3 dB 
MHz^^ cm^^ and y ~ 0.9 [13]. These values were assigned 
to the the annular region occupied by the acrylic in the 2D 
SOS maps co(r), density map po{v) and attenuation coefficient 
Q;o(r), respectively. The SOS value 1480 m s^^ and density 
value 1000 kg m"'^ of water were assigned elsewhere. Since 
we neglected the relatively weak acoustic attenuation due to 
the water bath and agar, Q!o(r) was also set to zero elsewhere. 

The experimental data were acquired from a cylindrically 
focused ultrasound transducer that had a central frequency 
of 2.25 MHz with a bandwidth of 70% [47]. The transducer 
was scanned along a circular trajectory of radius 95 mm, and 
20,000 temporal samples were measured at each transducer 
location at a sampling rate of 20 MHz. More details about 
the data acquisition can be found in Ref. [13]. In this study, 
images were reconstructed by use of PA signals recorded at 
200, 100 (over 180 degrees), and 50 transducer locations, 
which correspond to the full-view, limited-view, and few- 
view scanning geometry, respectively. The TR reconstruction 
method employed 20,000 temporal samples at each transducer 
location, while the iterative method employed 2,000 samples. 
The reference images were also reconstructed by use of the 
data obtained at 200 transducer locations when the acrylic 
cylinder was absent. Since the pencil lead phantom is expected 
to generate quasi-cylindrical waves and the morphology of the 
acoustic heterogeneity (the acrylic shell) was a cylinder, the 
cylindrical wave propagation can be approximated by the 2D 
PA wave equation. Accordingly, we employed a 2D imaging 
model in the experimental study, and all the reconstructions 
were performed on a grid of 512 x 512 pixels with a pitch of 
0.5 mm. The effects of shear wave propagation in the acrylic 
cylinder were neglected, which we expected to be of second- 
order importance compared to wavefield perturbations that 
arise from inhomogeneties in the SOS and density distributions 
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[48]. 



V. Simulation and experimental results 



A. Computer-simulations corresponding to different scanning 
geometries 

The reconstructed images corresponding to the three scan- 
ning geometries are displayed in Figs. 4 - 7. In each figure, the 
results in the top row correspond to use of the TR reconstruc- 
tion method, while the bottom row shows the corresponding 
results obtained by use of the iterative method. The profiles 
shown in each figure are along the 'Y'-axis indicated by the 
arrow in Fig. 4-(a). The red solid lines and blue dashed lines 
correspond to profiles through the phantom and reconstructed 
images, respectively. With the full-view scanning geometry, 
the TR method and the iterative method both produce accurate 
reconstructed images. However, with the few-view and the 
limited-view scanning geometries, the images reconstructed 
from the iterative method contain fewer artifacts and less 
noise than the TR results Also, the values of the images 
reconstructed from the iterative method are much closer to the 
values of the phantom than those produced by the TR method. 
The root mean square error (RMSE) between the phantom 
and the reconstructed images were also computed. The RMSE 
of images reconstructed by use of the TR method and the 
iterative method corresponding to noisy pressure data with the 
full-view, few-view, and limited-view scanning geometries are 
0.011, 0.042, 0.081 and 0.003, 0.007, 0.008, respectively. The 
computational time of the TR method was 1 .7 minutes, while 
the iterative method took approximately 10 minutes to finish 
20 iterations. 

B. Simulation results with errors in SOS and density maps 

Figure 8 shows the images reconstructed from noisy pres- 
sure data corresponding to the low contrast disc phantom in 
the case where SOS and density maps have no error The 
results corresponding to TR and iterative image reconstruction 
algorithms are shown in the top and bottom row, respectively. 
The RMSE corresponding to the time-reversal and the iterative 
results are 0.026 and 0.007, respectively. These results suggest 
that the iterative algorithm can more effectively reduce the 
noise level in the reconstructed images than the time-reversal 
algorithm. 

The images reconstructed by use of the SOS and density 
maps with errors are shown in Fig. 9. The image produced 
by the iterative method has cleaner background than the TR 
result, and the RMSE corresponding to the TR and the iterative 
results are 0.086 and 0.034, respectively. The boundaries of the 
disc phantoms also appear sharper in the image reconstructed 
by the iterative method as compared to the TR result. This can 
be attributed to the TV regularization employed in the iterative 
method. These results suggest that appropriately regularized 
iterative reconstruction methods can be more robust to the 
errors in the SOS and density maps than the TR method. 

^With the limited view scanning geometiy, we also implemented the iterated 
TR method [49], which produced images with fewer artifacts than the ordinary 
TR results, but the background was still not as clean as the iterative results. 
Given the limited space, those results were not included in this article. 
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Fig. 4. (a) and (c) are reconstructed images from noiseless data with full-view 
scanning geometry by use of the TR method and iterative method, respectively, 
(b) and (d) are the corresponding profiles along the 'Y'-axis indicated in panel 
(a). 
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Fig. 5. (a) and (c) are reconstructed images from the noisy pressure data with 
3% AWGN corresponding to the full-view scanning geometry by use of the TR 
method and iterative method, respectively, (b) and (d) are the corresponding 
profiles. 



C. 3D simulation results 

The 3D blood vessel phantom and the reconstructed images 
were visualized by the maximum intensity projection (MIP) 
method. Figure lO-(a) shows the phantom image, and Fig. 
lO-(b) and (c) display the images reconstructed by use of 
the TR method and the iterative method, respectively. They 
are all displayed in the same grey scale window. The RMSE 
corresponding to the TR and the iterative results are 0.018 
and 0.003, respectively. These results suggest that the iterative 
method is robust to the data incompleteness and the noise in 
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Fig. 6. (a) and (c) are reconstructed images from tlie noisy pressure data 
witli 3% AWGN corresponding to the few-view scanning geometry by use 
of tlie TR metliod and iterative method, respectively, (b) and (d) are the 
corresponding profiles. 



Fig. 8. (a) and (c) are reconstructed images with actual SOS and density 
maps by use of the TR method and iterative method, respectively, (b) and (d) 
are the con'esponding profiles along the 'Y'-axis indicated in panel (a). 
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Fig. 7. (a) and (c) are reconstructed images from the noisy pressure data 
with 3% AWGN corresponding to the limited-view scanning geometry by 
use of the TR method and iterative method, respectively, (b) and (d) are the 
corresponding profiles. 




(c) 



■ 0.5 



-^0 -20 -10 10 20 
Position Y (mm) 

(d) 



30 



Fig. 9. (a) and (c) are reconstructed images with SOS and density maps 
with errors by use of the TR method and iterative method, respectively, (b) 
and (d) are the corresponding profiles along the 'Y'-axis indicated in panel 
(a). 



the pressure data. The computational time of the TR method 
was approximately 6 minutes, while the iterative method with 
10 iterations required 110 minutes. 

D. Experimental results 

The images reconstructed from the experimental data are 
shown in Figs. 11 - 14. Figure 11 shows the image recon- 
structed with the full-view scanning geometry by use of the 
TR method (top row) and the iterative method (bottom row). 
Figure 11 -(a) and (c) display the reference images produced 



by each of the methods when the acrylic shell was absent. 
Figure 11 -(b) and (e) show the reconstructed images for the 
case when the acrylic shell was present. The RMSE between 
Fig. ll-(b), (d) and the reference images ll-(a), (c) are 0.003 
and 0.002, respectively. Figure 12-(a) and (c) show the images 
reconstructed with the few-view scanning geometry when the 
acrylic shell was present. The corresponding image profiles 
are displayed in Figure 12-(b) and (d). The profiles of Fig. 12- 
(a) and (c) along the 'Y'-axis were shown in Fig. 13, which 
shows that the iterative method produced higher resolution 
images than the TR method. This can be attritubed to the 
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Fig. 10. Maximum intensity projection renderings of the 3D pliantom (a), and the reconstructed 3D images by use of the TR method (b) and the iterative 
method (c). 
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Fig. 11. (a) and (b) are reconstructed images by use of the TR method 
from 200 views with acrylic shell absent and present, respectively, (c) and (d) 
are reconstructed images by use of the iterative method from 200 views with 
acrylic shell absent and present, respectively. 



Fig. 12. (a) and (c) are reconstructed images with data from 50 view angles 
over 360 degrees (acrylic shell present) by use of the TR method and iterative 
method, respectively, (b) and (d) are their corresponding profiles (dashed blue 
lines), where red solid lines are the profiles of the reference images in Fig. 
11 (a) and (c). 



TV regularization that mitigates model errors that arise, for 
example, by neglecting the shear wave and finite transducer 
aperture effects. The RMSE between Fig. 12-(b), (d) and their 
reference images are 0.005 and 0.002, , respectively. Figure 
14 displays the images reconstructed with the limited-view 
scanning geometry when the acrylic shell was present. The 
RMSE between Fig. 14-(a), (c) and their reference images 
are 0.007 and 0.003, respectively. These results show that 
the iterative algorithm can effectively compensate for the 
acoustic attenuation and mitigate artifacts and distortions due 
to incomplete measurement data. 

VI. Conclusion and discussion 

We proposed and investigated a full-wave approach to iter- 
ative image reconstruction in PACT with acoustically inhomo- 
geneous lossy media. An explicit formulation of the discrete 
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Fig. 13. The profiles of the reconstructed images in Fig. 12 along the 'Y'-axis 
indicated in Fig. 12(a). 



imaging model based on the k-space pseudospectral method 
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Fig. 14. (a) and (c) ai'e reconstructed images with data from 100 view angles 
over 180 degrees (acrylic shell present) by use of the TR method and iterative 
method, respectively, (b) and (d) are their con'esponding profiles (dashed blue 
lines), where red solid lines are the profiles of the reference images in Fig. 
11 (a) and (c). 



was described, and the details of implementing the forward 
and backprojection operators were provided. The matched 
operator pair was employed in an iterative image reconstruc- 
tion algorithm that sought to minimize a TV-regularized PLS 
cost function. The developed reconstruction methodology was 
investigated by use of both computer-simulated and experi- 
mental PACT measurement data, and the results demonstrated 
that the reconstruction methodology can effectively mitigate 
image artifacts due to data incompleteness, noise, finite sam- 
pling, and modeling errors. This suggests that the proposed 
image reconstruction method has the potential to be adopted 
in preclinical and clinical PACT applications. 

There remain several important topics to further investigate 
and validate the proposed iterative reconstruction method. It 
has been shown [20], [23] that the performance of recon- 
struction methods can be degraded when the SOS distribution 
satisfies a trapping condition [20], [23]. Therefore, future 
studies may include the investigation of numerical properties 
of the proposed image reconstruction method for cases in 
which the SOS distribution satisfies the trapping condition. 
Also, because the signal detectability is affected by the noise 
properties of an image reconstruction method, investigation 
of statistical properties of the iterative image reconstruction 
method is another important topic for future studies. Moreover, 
the proposed image reconstruction method can be further 
validated through additional experimental studies, and the 
quality of the produced images will be assessed by use of 
objective and quantitative measures. 
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APPENDix-A: Modeling transducer impulse 

RESPONSES 

An important feature of the proposed discrete PACT imag- 
ing model is that the transducer's impulse responses, including 
the spatial impulse response (SIR) and the acousto-electrical 
impulse response (EIR), can be readily incorporated into the 
system matrix. 

The SIR accounts for the averaging effect over the trans- 
ducer surface [50]-[52], which can be described as 



sW) ' 



(45) 



where ^^^^{rf ,mAt) is the averaged pressure at time t = 
mAt over the surface of the l-th transducer, S{rf) is the 
surface area of the l-th transducer centered at rf. 

In order to incorporate the SIR into the system matrix, we 
can divide the transducer surface into K small patches with 
equal area AS that is much less than the acoustic wavelength, 
so the integral in Eqn. 45 can be approximated by summation 
as 

AC 

>(rf,mAi)^;^, (46) 



/"^(rf,mAi) 



K 

fe=i 



or in the equivalent matrix form 



(47) 



where r^*^ denotes the center of the fc-th patch of the ^-th 



transducer, AS* is the patch area, 7 



SIR 



AS 

s(rf: 



1) is a 



1 X K vector, pj„ = {p{^h '^tAi), • • • , p{rf , rjiAt))^ denotes 
the acoustic pressure at patches of l-th transducer at time 
mAt. Here for simplicity, we assume all the transducers are 
divided into K patches with equal area AS', and it is readily 
to extend to general cases where l-th transducer is divided into 
Ki patches with area of ASik. 

Recalling the measured pressure data and p defined for 
point-Uke transducer, we can redefine Pm as a KL x 1 vector 
that represents the acoustic pressure at patches of transducers 
with finite area at time t = mAt as 



(48) 



The corresponding p can be redefine as a KLM x 1 vector 
denoting the measured pressure data corresponding to all 
transducer and temporal samples as 



Po 



Pm-1 



(49) 



The averaged pressure measured by all transducer and 
temporal samples can be defined as the LM x 1 vector 



-SIR _ 



-SIR 

Pm-i. 



(50) 
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where the L x 1 vector 



-SIR 



According to Eqn. 47, p and p^"^ can be related as 



^SIR 

where the KLM x LM matrix 
„SIR 



pSIR^ 



-,SIR _ 
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SIR 



(51) 



(52) 



(53) 



The EIR models the electrical response of the piezoelectric 
transducer With the assumption that the transducer is a linear 
shift invariant system with respect to the input averaged pres- 
sure time sequence, the output voltage signal is the convolution 
result of the input and the EIR. 

For simplicity, the transducers are assumed to process 
identical EIR, and let h"^ = (/if,-- - jhj)"^ be the discrete 
samples of the EIR. The input averaged pressure time se- 
quence of the l-th transducer can be defined as a L x 1 
vector p^iR = (p^"^(rf , 0), - - - ,p^"^(rf , (M - 1)A<))T. Then 
the output voltage signal pf- of the Z-th transducer can be 
expressed as a (J + i\/— 1) x 1 vector 



i^" * PSIR, 



(54) 



where * denotes discrete linear convolution operation, which 
can be constructed as a matrix multiplication by converting 
one of the operands into the corresponding Toeplitz matrix. 

The output voltage signals of all transducers p''* = 
(Pi^i ■ ■ ■ J P^)^ can then be computed as 



-IR 



pEIR^SIR 



where the L{J + M — 1) x LAI matrix 

„EIR1 

-,EIR 
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(55) 
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and 7 

as 



EIR 



is a (J + ii/ — 1) X LM Toeplitz-like matrix defined 
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(57) 

By use of Eqns. (36), (52), and (55), it is readily found that 

,IR _ T-t^lKT-^lKcrr,, .TlToPo, (58) 



■nEIRT-iSIRorxi 



and the corresponding system matrix that incorporates the 
transducer impulse responses is found to be 

jjiR ^ r^iRp^iRsTM-i • • • TiTo. (59) 



APPENDix-B: Implementation of theFISTA 

ALGORITHM FOR PACT 

Equation (43) was solved iteratively whose pseudocodes are 
provided in Alg. 1, where 'Lip' is the Lipschitz constant of 
the operator 2H'^H [42]. 

Algorithm 1 Solver of the optimization problem defined by 
Eqn. (43) 

Input: p, Pq°\ a. Lip 
Output: po 

1: ^ 1; (T^^^ ^ p[)°^ {Set the initial guess (The 
zero initial guess was employed in all the studies in this 
article)} 
2: for C = 1 to Z do 



'^0 



P),2A/Lip) 



Lip 



JC+i) 
end for 

(Z) 

Po ^ Po 



Note that we extended the FISTA algorithm described in 
Ref. [42] to 3D. The function 'F_Dnoise' in Alg. LLine 3 
solves a de-noising problem defined as: 



X = argmin||y — x| 

x>0 



-/3|x|tv, 



where /3 = 2 A/Lip and 



2 
Lip 



y = p-7^HT(HrT|,^)-p). 



(60) 



(61) 



It has been demonstrated that Eqn. (60) can be solved effi- 
ciently [42], and the pseudocodes are provided in Alg. 2. 

Algorithm 2 Solver of the de-noising problem defined by Eqn. 

(60) 

Input: y, /3 
Output: X 

a(i),b(i),c(i)] ^ 

O(Afi-l) X xAf.-i ) X (JV2-I) xiVa ' X JVa x(Af3-l)] 

d(0),e(0),f(0)] ^ 

O(Afi-l) X Afa xWa J X (JV2-I) xWs 7 OaTj X JV2 X (W.,-!)] 
= 1 

2: for C = 1 to Z do 

3: [d(0,e(C),f(C)] ^ Pp{[a(';),b(«,c«)] + 

4: ^ 1 + o.5vTT4(FO)2 

5: [a(C+i),b(C+i),c(«+i)] ^ - l)/t(C+i)[d(0 - 

d(f-i),e(C)-e(^-i),f(«)-f('^-i)] 
6: end for 

7: x^7'c{y-A7'({d(^),e(^),f(^)}} 
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The four operators Vi Vc, Vf and Vp in Alg. 2 are defined 
as follows: 

Vl : R(^i-1)x^2X7V3 x]]jJViX(Af2-l)xW3 ^jjATiXJVaxCAfa-l) _^ 
^NixN2y.N3 ^ 

[7'i{a,b,c}] 

L ' ' J J m ,712 ,^l3 

[^]ni,n2,n3 H" [^]ni,n2,n3 \S^Ti\ ,712 ,1^-^ 

[^]ni — 1 ,712 ,713 ['^] ^T'l i^T'^ — 1 i^T'S ~^ [^] ni .712 , 713 — 1 

for ni = I,-- - ,7Vi, 712 = I,-- - ,^2, "3 = I,-- - ,^^3, 

(62) 

where we assume [a]o,„2,„3 = [a]Ari,ri2,n3 = [b],ii,o,n3 = 

[b],ii,Ar2,n3 = [c]rai,n2,0 = [c]ni ,n2 ,-/V3 = 0. 
Vc ■ ^'^^ XN2XN3 _^ j^Wi XN2XN3 

['Pc{^}]m^n2,n3 = max{0, [x]„,,„,,„3 .} (63) 

-pT . ^NixN2XN3 _j. Ju(Wi-l)xJV2XJV3 ^ ]^Afix(Ar2-l)xAf3 ^ 

]]^AfixAf2x(Af3-i) ,^^g denote the input and output matrices 
by y and (a, b, c) respectively, we have 

,ri2 ,Ti3 = [y]ni,n2,"3 [y] )ii +I,n2 ,113 ) 

for ni = I,-- - ,A^i - l,7i2 = I,-- - ,^2,"-3 = I,-- - ,^3 

[t)]ni,n2,n3 — [y]ni ,712,713 [y]ni ,712 + 1,713 ; 

for rii = I,-- - ,A^i,7T,2 = I,-- - ,iV2- l,?i3 = l,-" ,^3 

[c]ni,7l2,7l3 = [y]ni,7l2,n3 [y]7ll ,712,713 + 1 1 

for ni = 1, • • • , A^i,ri2 = 1, • • • ,iV2,?^3 = 1, • • • ,^3 - 1- 

(64) 

•P ■.Ri^^-^)'^^2XN3 y,-^NiX{N2-l)xN3 y,^NiXN2X{N3-l) _^ 

^iNi-l)xN2XN3 y.-^Nix{N2-l)xN3 ^ R^i x ^2 X (Ar3-1) ^ If We 

denote the input and output matrices by (a, b, c) and (d, e, f ) 
respectively, we have 
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[t>] 111, 712, 713 
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' Y^t'^-l^^l -^2 ,713 ^ [b]^^ ^^2,7l3 ^ 
[^]ni ,712 ,713 


L"^J71i ,712.713 , 


^[17] 
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\cP } 

L^j7li ,712 ,713 J 


[18] 



(65) 



where ni = 1, • • • , iVi, 712 = 1, • • • , iV2, "-3 = 1, ' ' ' , ^3. 

and we assume [a]o,,i2,7i3 = [a]Afi,7i2,7i3 = [b]7ii,o,7i3 

[b]7ii ,Af2,n3 = [c]ni, 712,0 = ['^]7ll ,712 ,-/V3 = 0. 
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